source('00_util_scripts/mod_bulk.R')
source('00_util_scripts/mod_bplot.R')
library(readxl)
library(clusterProfiler)

proj.nm <- 'mission/FPP/'

# 293 ----------
hek293 <- read_csv('mission/FPP/thapsigargin_upr/GSE293666_HEK293_IRE1_activators_raw_counts_all_samples.csv.gz')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

fct_mva <- tibble(gene = kegg_mva, order = fct_inorder(gene) |> fct_rev())

mva293 <- hek293 |>
  pivot_longer(-1) |>
  calc_tpm(sample = name, abundance = value) |>
  select(name, Name, tpm) |>
  filter(Name %in% kegg_mva, str_detect(name, 'DMSO|Tg'),
         tpm > 0)

mva293 |>
  left_join(fct_mva, join_by(Name == gene)) |>
  tidyplot(name, order, color = tpm) |>
  add_heatmap(scale = 'row') +
  labs(fill = 'Z-score', x = 'Sample', y = 'Gene',
       subtitle = 'GSE293666',
       title = 'MVA pathway HEK293 after 6h treatment of 0.5uM TG')

publish_source_plot('293_tg_MVA')

tdb293 <- hek293 |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, 'Tg|DMSO')) |>
  na.omit() |>
  tidybulk(.sample = name, .transcript = Name,
           .abundance = value) |>
  quick_process_bulk(group = group, skip_scale = T)

res293 <- tdb293 |>
  test_differential_abundance(~0+group, omit_contrast_in_colnames = T,
                              contrasts = 'groupTg-groupDMSO') |>
  pivot_transcript()

res293 |>
  filter(Name %in% kegg_mva) |>
  ggplot(aes(logFC, Name, fill = FDR < .05)) +
  geom_col() +
  scale_fill_manual(values = c('grey','red'))

res293 |>
  write_csv('mission/FPP/thapsigargin_upr/293hek_tg_vs_dmso_deg.csv')

gsego293 <- res293 |>
  filter(FDR < .05) |>
  pull(logFC, name = Name) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db',keyType = 'SYMBOL',
        pvalueCutoff = 1)

tg293_upr_nes <- gsego293@result |>
  as_tibble() |>
  filter(str_detect(Description, 'unfold')) |>
  select(ID, Description, NES, qvalue, leading_edge)

# skin fibroblast -------
tg_fibro <- read_tsv('mission/FPP/thapsigargin_upr/GSE27349.top.table.tsv')

geo_fibro <- pluck_geo('GSE27349', annot_gpl = T)

fData(geo_fibro) |> head()

pData(geo_fibro) |> head()

acc_mva <-
bitr(kegg_mva, fromType = 'SYMBOL', toType = 'REFSEQ', OrgDb = 'org.Hs.eg.db')

tg_fibro |>
  inner_join(acc_mva, join_by(GB_ACC == REFSEQ))

symb_fibro <-
tg_fibro$GB_ACC |>
  bitr(toType = 'SYMBOL', fromType = 'ACCNUM', OrgDb = 'org.Hs.eg.db')

symb_fibro |>
  as_tibble() |>
  filter(SYMBOL %in% kegg_mva)

# hepalike hESC ------
hesc <- 
  read_excel('mission/FPP/thapsigargin_upr/GSE284753_CP20M_240620_A00154_1531_AHJ273DSXC_LTR.xlsx', .name_repair = 'universal')

colnames(hesc) <- colnames(hesc) |>
  str_remove('\\.R\\d.+')

hesc_wt <- hesc |>
  select(Gene.Name, LTR1, LTR2, LTR5, LTR6, LTR9, LTR10)

hesc_tdb <- hesc_wt |>
  pivot_longer(-1) |>
  mutate(suffix = str_extract(name, '\\d+') |> as.integer(),
         group = ifelse(suffix %% 2 == 0, 'Thaps', 'DMSO')) |>
  tidybulk(.sample = name,
           .transcript = Gene.Name,
           .abundance = value)

hesc_tdb <- hesc_tdb |>
  quick_process_bulk(group = group, skip_scale = T)

hesc_res <- hesc_tdb |>
  test_differential_abundance(~ 0 + group, omit_contrast_in_colnames = T,
                              contrasts = 'groupThaps-groupDMSO')

hesc_res |>
  pivot_transcript() |>
  filter(Gene.Name %in% kegg_mva) |>
  ggplot(aes(logFC, Gene.Name, fill = FDR < .05)) +
  geom_col() +
  theme_bw() +
  scale_fill_manual(values = c('grey','red'), labels = c('NS', 'Upregulated')) +
  labs(fill = 'Significance',
       subtitle = 'Traini, et al. Stem Cell Res Ther. 2025',
       title = 'MVA pathway in H1 human ESC treated with 8h treatment of 1uM TG')

hesc_tdb |>
  calc_tpm(group = group, abundance = value) |>
  mutate(sample = str_c(group, suffix),
         Gene.Name, tpm, .keep = 'none') |>
  filter(Gene.Name %in% kegg_mva) |>
  as_tibble() |>
  heatmap(Gene.Name, sample, tpm, scale = 'row',
          palette_value = c("#67A9CF", "#F7F7F7", "#EF8A62")) +
  labs(x = 'Sample')

hesc_tdb |>
  calc_tpm(sample = name, abundance = value) |>
  mutate(sample = str_c(group, suffix),
         Gene.Name, tpm, .keep = 'none') |>
  filter(Gene.Name %in% kegg_mva) |>
  left_join(fct_mva, join_by(Gene.Name == gene)) |>
  tidyplot(sample, order, color = tpm) |>
  add_heatmap(scale = 'row') +
  labs(fill = 'Z-score',
       subtitle = 'Traini, et al. Stem Cell Res Ther. 2025',
       title = 'MVA pathway in H1 human ESC treated with 8h treatment of 1uM TG')

publish_source_plot('hESC_tg_MVA')

hesc_res <- hesc_res |>
  pivot_transcript() |>
  select(Gene.Name, logFC, PValue, FDR) |>
  write_csv('mission/FPP/thapsigargin_upr/hesc_tg_vs_dmso_deg.csv')

hesc_gsego <- hesc_res |>
  filter(FDR < .05) |>
  pull(logFC, name = Gene.Name) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db',keyType = 'SYMBOL',
        pvalueCutoff = 1)

tg_hesc_upr_nes <- hesc_gsego@result |>
  as_tibble() |>
  filter(str_detect(Description, 'unfold')) |>
  select(ID, Description, NES, qvalue, leading_edge)

tg_hesc_upr_nes  

## ER UPR NES sum-up ----------
list('TG_hESC' = tg_hesc_upr_nes, 'TG_293HEK' = tg293_upr_nes,
     'SA_mouse_KC' = sa_kc_upr_nes, 'UVB_skin' = uvb_skin_upr_nes) |>
  list_rbind(names_to = 'condition') |>
  filter(ID == 'GO:0030968') |>
  mutate(condition = fct_reorder(condition, NES)) |>
  ggplot(aes(NES, condition, fill = qvalue)) +
  geom_col() +
  theme_pubr(legend = 'right') +
  scale_fill_distiller(palette = 'Reds') +
  labs(title = 'Enrichment of ER UPR pathway gene in UPR inducing conditions',
       x = 'Normalized enrichment score', fill = 'FDR')

### chaperones logfc ---------
chaperon <- read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

sa_kc_chap <- kera_v3_deg |>
  filter(str_to_upper(gene) %in% chaperon$SYMBOL) |>
  mutate(FDR = p_val_adj, logFC = avg_log2FC, gene = str_to_upper(gene),
         Condition = 'SA_mouse_KC', .keep = 'none')

uvb_skin_chap <- res.sk6 |>
  filter(Symbol %in% chaperon$SYMBOL) |>
  mutate(gene = Symbol, Condition = 'UVB_skin', FDR, logFC, .keep = 'none')

tg_hesc_chap <- hesc_res |>
  filter(Gene.Name %in% chaperon$SYMBOL) |>
  mutate(gene = Gene.Name, Condition = 'TG_hESC', FDR, logFC, .keep = 'none')

tg_293_chap <- res293 |>
  filter(Name %in% chaperon$SYMBOL) |>
  mutate(gene = Name, Condition = 'TG_293HEK', FDR, logFC, .keep = 'none')

sum_chap <-
  list(sa_kc_chap, uvb_skin_chap, tg_hesc_chap, tg_293_chap) |>
  list_rbind() 

sum_chap |>
  tidyplot(gene, Condition, color = logFC, width = 100) |>
  add_heatmap() |>
  adjust_colors(new_colors = colors_diverging_blue2red,
                values = pretty_distiller(sum_chap$logFC)) |>
  add_title('UPR-associated chaperone logFC in UPR inducing conditions')

proj.nm <- 'mission/FPP/'

publish_source_plot('SA.TG.UVB.UPR.chaperone.heatmap', width = 100)

sum_chap |>
  mutate(features.plot = gene, id = Condition, avg_log2fc = logFC,
         p_val_adj = FDR) |>
  BubblePlot()

sum_chap |>
  filter(gene %in% c('CALR','HSP90B1','PARK7','SGTA','HSPA5')) |>
  ggplot(aes(gene, logFC, fill = FDR)) +
  geom_col() +
  facet_wrap(~Condition, scales = 'free_y') +
  theme_bw() +
  scale_fill_distiller(palette = 'Reds') +
  labs(title = 'UPR-associated chaperone logFC in UPR inducing conditions')

sum_chap |>
  filter(gene %in% c('CALR','HSP90B1','PARK7','SGTA','HSPA5')) |>
  ggplot(aes(Condition, logFC, fill = gene)) +
  geom_col(position = 'dodge') +
  scale_fill_viridis_d(option = 'turbo', begin = .1) +
  theme_bw() +
  labs(title = 'UPR-associated chaperone logFC in UPR inducing conditions')

### chaperone violin/barplot ------------
upr_chap5 <- c('CALR','HSP90B1','PARK7','SGTA','HSPA5')

tdb293 |>
  calc_tpm(sample = name, abundance = value) |>
  mutate(log_tpm = log1p(tpm)) |>
  filter(Name %in% upr_chap5) |>
  ggplot(aes(group, log_tpm, fill = group)) +
  stat_summary(geom = 'col', fun = 'mean') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_normal', width = .5) +
  facet_wrap(vars(Name), nrow = 1, scales = 'free_y') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_hue(direction = -1) +
  labs(title = 'Tg treated 293HEK cells', y = 'log(tpm+1)') +
  theme_jpub

publish_source_plot('TG.293.upr.5chaperone.barplot', width = 100)

tdb293 <- read_csv('mission/FPP/results/TG.293.upr.5chaperone.barplot.csv')

tdb293 |>
  filter(group == 'DMSO') |>
  summarise(ref_ltpm = log1p(mean(tpm)), .by = Name) |>
  right_join(tdb293) |>
  mutate(Foldchange = log_tpm/ref_ltpm) |>
  ggplot(aes(group, Foldchange, fill = group)) +
  stat_summary(geom = 'col', fun = 'mean') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_normal', width = .5) +
  facet_wrap(vars(Name), nrow = 1) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_hue(direction = -1) +
  labs(title = 'Tg treated 293HEK cells', y = 'Fold change') +
  theme_jpub

publish_pdf('TG.293.upr.5chaperone.barplot.fc.pdf', width = 100)

hesc_tdb |>
  calc_tpm(sample = name, abundance = value) |>
  mutate(log_tpm = log1p(tpm), group = ifelse(group == 'DMSO', 'DMSO', 'Tg')) |>
  filter(Gene.Name %in% upr_chap5) |>
  ggplot(aes(group, log_tpm, fill = group)) +
  stat_summary(geom = 'col', fun = 'mean') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_normal', width = .5) +
  facet_wrap(vars(Gene.Name), nrow = 1, scales = 'free_y') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_hue(direction = -1) +
  labs(title = 'Tg treated hESC cells', y = 'log(tpm+1)') +
  theme_jpub

publish_source_plot('TG.hESC.upr.5chaperone.barplot', width = 100)

tdb.sk <- read_csv('mission/FPP/uvb/GSE148535.human.skin.bulk.csv') |>
  tidybulk(.sample = title,
           .transcript = Symbol,
           .abundance = value)

tdb.sk |>
  filter(group != '24_hour') |>
  mutate(log_tpm = log1p(tpm)) |>
  filter(Symbol %in% upr_chap5) |>
  ggplot(aes(group, log_tpm, fill = group)) +
  stat_summary(geom = 'col', fun = 'mean') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_normal', width = .5) +
  facet_wrap(vars(Symbol), nrow = 1, scales = 'free_y') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  scale_fill_hue(direction = -1) +
  labs(title = 'UVB treated skin tissue', y = 'log(tpm+1)') +
  theme_jpub

publish_source_plot('UVB.skin.upr.5chaperone.barplot', width = 100)

## SREBF2 target genes ------------
srebf2_hs_acti <-
read_tsv('00_util_scripts/ref/trrust_rawdata.human.tsv',
         col_names = c('TF','target','type','ref')) |>
  filter(TF == 'SREBF2', type != 'Repression')

srebf2_mm_acti <-
  read_tsv('00_util_scripts/ref/trrust_rawdata.mouse.tsv',
           col_names = c('TF','target','type','ref')) |>
  filter(TF == 'Srebf2', type != 'Repression')

srebf2_mm_acti

sa_kc_deg <- kera_v3_deg |>
  mutate(FDR = p_val_adj, logFC = avg_log2FC, gene = str_to_upper(gene),
         Condition = 'SA_mouse_KC', .keep = 'none')

uvb_skin_deg <- res.sk6 |>
  mutate(gene = Symbol, Condition = 'UVB_skin', FDR, logFC, .keep = 'none')

tg_hesc_deg <- hesc_res |>
  mutate(gene = Gene.Name, Condition = 'TG_hESC', FDR, logFC, .keep = 'none')

tg_293_deg <- res293 |>
  mutate(gene = Name, Condition = 'TG_293HEK', FDR, logFC, .keep = 'none')

sum_upr_deg <-
  list(sa_kc_deg, uvb_skin_deg, tg_hesc_deg, tg_293_deg) |>
  list_rbind()

sum_upr_srebf2 <- sum_upr_deg |>
  filter(gene %in% srebf2_hs_acti$target)

sum_upr_srebf2 |>
  tidyplot(gene, Condition, color = logFC, width = 50) |>
  add_heatmap() |>
  adjust_colors(new_colors = colors_diverging_blue2red,
                values = pretty_distiller(sum_chap$logFC)) |>
  add_title('SREBF2 target genes logFC in UPR inducing conditions')

term2gene_srebf2 <- srebf2_hs_acti |>
  mutate(term = 'SREBF2 targets', gene = target, .keep = 'none') |>
  distinct(gene, .keep_all = T)

sa_kc_srebf2_nes <- sa_kc_deg |>
  #filter(FDR<.05) |>
  pull(logFC, name = gene) |>
  sort(decreasing = T) |>
  GSEA(TERM2GENE = term2gene_srebf2, pvalueCutoff = 1, pAdjustMethod = 'none',
       minGSSize = 3)

sa_kc_srebf2_nes |>
  gseaplot2('SREBF2 targets')

sa_kc_deg |>
  #filter(FDR < .05) |>
  mutate(rank = rank(-logFC)) |>
  filter(gene %in% kegg_mva)

# prepare TRAPT inputs ---------
## SA KC --------
sa_kc_deg |>
  filter(FDR < .05, logFC > 1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/sa_kc_upfc0.txt', col_names = F)

sa_kc_trapt <-
  read_delim('~/append-ssd/learn/learn_trapt/sa_kc_upfc0_out/TR_detail.txt')

sa_kc_trapt |>
  filter(tr_base == 'SREBF2') |>
  ggplot(aes(Source, `TR activity`)) + geom_boxplot() +
  labs(title = 'SA KC SREBF2')

sa_kc_deg |>
  filter(FDR < .05, logFC < -1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/sa_kc_downfc1.txt', col_names = F)

sa_kc_trapt_down <-
  read_delim('~/append-ssd/learn/learn_trapt/sa_kc_downfc1_out/TR_detail.txt')

list('SA'=sa_kc_trapt, 'PBS'=sa_kc_trapt_down) |>
  list_rbind(names_to = 'group') |>
  filter(tr_base == 'SREBF1', Source == 'Cistrome') |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  labs(title = 'SA KC SREBF2')

## uvb skin ------------
uvb_skin_deg |>
  filter(FDR < .05, logFC > 1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/uvb_skin_upfc1.txt', col_names = F)

uvb_skin_trapt <-
  read_delim('~/append-ssd/learn/learn_trapt/uvb_skin_upfc1_out/TR_detail.txt')

uvb_skin_trapt |>
  filter(tr_base == 'SREBF2') |>
  ggplot(aes(Source, `TR activity`)) + geom_boxplot() +
  labs(title = 'UVB skin SREBF2')

uvb_skin_deg |>
  filter(FDR < .05, logFC < -1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/uvb_skin_downfc1.txt', col_names = F)

uvb_skin_trapt_down <-
  read_delim('~/append-ssd/learn/learn_trapt/uvb_skin_downfc1_out/TR_detail.txt')

list('SA'=uvb_skin_trapt, 'PBS'=uvb_skin_trapt_down) |>
  list_rbind(names_to = 'group') |>
  filter(tr_base == 'SREBF1', Source == 'Cistrome') |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  labs(title = 'UVB skin SREBF1')

## TG 293 -------------
tg_293_deg |>
  filter(FDR < .05, logFC > 1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/tg_293_upfc1.txt', col_names = F)

tg_293_trapt <-
  read_delim('~/append-ssd/learn/learn_trapt/tg_293_upfc1_out/TR_detail.txt')

tg_293_trapt |>
  filter(tr_base == 'SREBF2') |>
  ggplot(aes(Source, `TR activity`)) + geom_boxplot() +
  labs(title = 'TG 293HEK SREBF2')

tg_293_deg |>
  filter(FDR < .05, logFC < -1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/tg_293_downfc1.txt', col_names = F)

tg_293_trapt_down <-
  read_delim('~/append-ssd/learn/learn_trapt/tg_293_downfc1_out/TR_detail.txt')

list('TG'=tg_293_trapt, 'DMSO'=tg_293_trapt_down) |>
  list_rbind(names_to = 'group') |>
  filter(tr_base == 'SREBF1', Source == 'Cistrome') |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  labs(title = 'TG 293HEK SREBF1')

## TG hESC ------------
tg_hesc_deg |>
  filter(FDR < .05, logFC > 1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/tg_hesc_upfc1.txt', col_names = F)

tg_hesc_trapt <-
  read_delim('~/append-ssd/learn/learn_trapt/tg_hesc_upfc1_tr_out/TR_detail.txt')

tg_hesc_trapt |>
  mutate(z_score = percent_rank(`TR activity`)) |>
  filter(tr_base == 'SREBF2') |>
  ggplot(aes(Source, z_score)) + geom_boxplot() +
  labs(title = 'TG hESC SREBF2')

tg_hesc_deg |>
  filter(FDR < .05, logFC < -1) |>
  select(gene) |>
  write_csv('~/append-ssd/learn/learn_trapt/tg_hesc_downfc1.txt', col_names = F)

tg_hesc_trapt_down <-
  read_delim('~/append-ssd/learn/learn_trapt/tg_hesc_downfc1_out/TR_detail.txt')

list('TG'=tg_hesc_trapt, 'DMSO'=tg_hesc_trapt_down) |>
  list_rbind(names_to = 'group') |>
  filter(tr_base == 'SREBF1', Source == 'Cistrome') |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  labs(title = 'TG hESC SREBF1')

## sum up ----------
trapt_sum <- 
  list('SA_KC'=sa_kc_trapt, 'UVB_skin'=uvb_skin_trapt,
     'TG_293HEK'=tg_293_trapt, 'TG_hESC'=tg_hesc_trapt) |>
  list_rbind(names_to = 'Condition') |>
  #filter(Source == 'ENCODE') |>
  mutate(rank_score = percent_rank(`TR activity`), .by = Condition) |>
  select(Condition, rank_score, tr_base, `TR activity`)

trapt_sum_down <- 
  list('SA_KC'=sa_kc_trapt_down, 'UVB_skin'=uvb_skin_trapt_down,
       'TG_293HEK'=tg_293_trapt_down, 'TG_hESC'=tg_hesc_trapt_down) |>
  list_rbind(names_to = 'Condition') |>
  #filter(Source == 'ENCODE') |>
  mutate(rank_score = percent_rank(`TR activity`), .by = Condition) |>
  select(Condition, rank_score, tr_base, `TR activity`)

list('UPR'=trapt_sum, 'Control'=trapt_sum_down) |>
  list_rbind(names_to = 'group') |>
  filter(str_detect(tr_base, 'SREBF1')) |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  facet_wrap(~Condition, nrow = 1) +
  labs(title = 'Predicted SREBF1 TF activity', subtitle = 'TRAPT Cistrome',
       y = 'TF activity percent rank')

list('UPR'=trapt_sum, 'Control'=trapt_sum_down) |>
  list_rbind(names_to = 'group') |>
  filter(str_detect(tr_base, 'SREBF2')) |>
  ggplot(aes(group, `TR activity`)) +
  geom_boxplot() +
  facet_wrap(~Condition, nrow = 1) +
  labs(title = 'Predicted SREBF2 TF activity', subtitle = 'TRAPT Cistrome',
       y = 'TF activity percent rank')

trapt_sum |>
  filter(tr_base == 'SREBF2') |>
  ggplot(aes(Condition, z_score)) +
  geom_boxplot() +
  labs(title = 'Predicted SREBF2 TF activity', subtitle = 'TRAPT Cistrome',
       y = 'TF activity percent rank') +
  theme_bw(base_size = 6, base_family = 'ArialMT')

publish_source_plot('trapt.cistrome.upr.condition.srebf2')

list('SA_KC'=sa_kc_trapt, 'UVB_skin'=uvb_skin_trapt,
     'TG_293HEK'=tg_293_trapt, 'TG_hESC'=tg_hesc_trapt) |>
  list_rbind(names_to = 'Condition') |>
  filter(Source == 'Cistrome') |>
  mutate(z_score = percent_rank(`TR activity`), .by = Condition) |>
  filter(tr_base == 'SREBF1') |>
  ggplot(aes(Condition, z_score)) +
  geom_boxplot() +
  labs(title = 'Predicted SREBF1 TF activity', subtitle = 'TRAPT Cistrome',
       y = 'TF activity percent rank') +
  theme_bw(base_size = 6, base_family = 'ArialMT')

publish_source_plot('trapt.Cistrome.upr.condition.srebf1')

list('SA_KC'=sa_kc_trapt, 'UVB_skin'=uvb_skin_trapt,
     'TG_293HEK'=tg_293_trapt, 'TG_hESC'=tg_hesc_trapt) |>
  list_rbind(names_to = 'Condition') |>
  filter(tr_base == 'SREBF2', Source == 'Cistrome') |>
  ggplot(aes(Condition, `TR activity`)) +
  stat_mean(geom = 'col') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_ci') +
  labs(title = 'Predicted SREBF2 TF activity', subtitle = 'TRAPT Cistrome') +
  theme_bw()

## random test TRAPT --------
ensembl_symbol_mapping |>
  filter(ref_genome == 'hg38') |>
  distinct(transcript) |>
  sample_n(1219) |>
  write_csv('~/append-ssd/learn/learn_trapt/random_gene1219.txt', col_names = F)

trapt_random <-
  read_delim('~/append-ssd/learn/learn_trapt/random_gene1219_out/TR_detail.txt')

trapt_random |>
  filter(Source == 'Cistrome', str_detect(tr_base, 'SREBF')) |>
  ggplot(aes(tr_base, `TR activity`)) +
  geom_boxplot()

trapt_random |>
  filter(Source == 'Cistrome') |>
  mutate(rank = rank(`TR activity`), .keep = 'used') |>
  ggplot(aes(rank, `TR activity`)) +
  geom_point()

# LISA ------------
lisa_upr <-
list.files('~/append-ssd/learn/learn_trapt/', 'basic', full.names = T) |>
  read_tsv(id = 'file')

## SREBF ----------
lisa_upr |>
  mutate(file = str_extract(file, '(?<=motif_).+(up|down)'),
         condition = str_remove(file, '_(up|down)'),
         type = str_extract(file, 'up|down')) |>
  filter(str_detect(factor, 'SREBF.$')) |>
  ggplot(aes(condition, -log10(summary_p_value), fill = type)) +
  geom_col(position = 'dodge') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_pubr() +
  facet_wrap(~factor) +
  scale_fill_manual(values = c('royalblue','tomato')) +
  labs(title = 'LISA predicted TF activity in UPR inducing condition',
       x = 'Condition')

## ATF4 ----------
lisa_upr |>
  mutate(file = str_extract(file, '(?<=motif_).+(up|down)'),
         condition = str_remove(file, '_(up|down)'),
         type = str_extract(file, 'up|down')) |>
  filter(str_detect(factor, '(ATF4|ATF6|XBP1)')) |>
  ggplot(aes(condition, -log10(summary_p_value), fill = type)) +
  geom_col(position = 'dodge') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_pubr() +
  facet_wrap(~factor) +
  scale_fill_manual(values = c('royalblue','tomato')) +
  labs(title = 'LISA predicted TF activity in UPR inducing condition',
       x = 'Condition')

# xcore ------------
library(xcore)
library(xcoredata)

s2f <- symbol2fantom()
remap_ref <- remap_promoters_f5()
chipat_ref <- chip_atlas_promoters_f5()

hek293 <- read_csv('mission/FPP/thapsigargin_upr/GSE293666_HEK293_IRE1_activators_raw_counts_all_samples.csv.gz')

tg293 <- hek293 |>
  select(matches('Name|DMSO|Tg')) |>
  column_to_rownames('Name') |>
  mutate(across(everything(), as.integer)) |>
  as.matrix()

hesc_int <- hesc_tdb |>
  tidybulk::aggregate_duplicates(.sample = name, .transcript = Gene.Name, .abundance = value) |>
  pivot_wider(names_from = name, values_from = value, id_cols = Gene.Name)
  
hesc_int <- hesc_int |>
  column_to_rownames('Gene.Name') |>
  mutate(across(everything(), as.integer)) |>
  as.matrix()

hesc_fntm <- hesc_int |> translateCounts(s2f)

hesc_design <- tibble(sample = colnames(hesc_int),
                      DMSO = str_ends(sample, '1|5|9'),
                      Tg = !DMSO) |>
  column_to_rownames('sample') |>
  as.matrix()

fantom293 <- tg293 |> translateCounts(s2f)

design293 <- tibble(sample = colnames(tg293),
                    DMSO = str_count(sample, 'DMSO'),
                    Tg = str_count(sample, 'Tg')) |>
  column_to_rownames('sample') |>
  as.matrix()

mae <- hesc_fntm |>
  prepareCountsForRegression(hesc_design, base_lvl = 'DMSO')

mae <- fantom293 |>
  prepareCountsForRegression(design = design293, base_lvl = "DMSO")

# it is ok to use a subset of sig ref
mae <- addSignatures(mae, remap = chipat_ref)

mae <- filterSignatures(mae, min = 0.05, max = 0.95)

# register parallel backend
# 16 worker cost ~2min for total remap ref
doParallel::registerDoParallel(16L)
BiocParallel::register(BiocParallel::DoparParam(), default = TRUE)

# take 5min+ to model chip-atlas in 16 p
res <- modelGeneExpression(
  mae = mae,
  xnames = "remap",
  nfolds = 5)

tg_293_xcore <- res$results$remap |>
  as_tibble()

tg_293_xcore_chipatlas <- res$results$remap |>
  as_tibble()

tg_293_xcore |>
  write_csv('mission/FPP/thapsigargin_upr/tg_293_xcore_TF.csv')

tg_hesc_xcore <- res$results$remap |>
  as_tibble() |>
  write_csv('mission/FPP/thapsigargin_upr/tg_hesc_xcore_TF.csv')

tg_293_xcore |>
  filter(str_detect(name, 'SREBF')) |>
  mutate(factor = str_extract(name, 'SREBF.')) |>
  summarise(z_score = max(z_score), .by = factor) |>
  ggplot(aes(factor, z_score)) +
  geom_col()

tg_hesc_xcore |>
  filter(str_detect(name, 'SREBF')) |>
  mutate(factor = str_extract(name, 'SREBF.')) |>
  summarise(z_score = max(z_score), .by = factor) |>
  ggplot(aes(factor, z_score)) +
  geom_col()

tg_293_xcore |>
  filter(str_detect(name, 'ATF4|ATF6|XBP1')) |>
  ggplot(aes(z_score, name)) +
  geom_col()

tg_293_xcore_chipatlas |>
  filter(str_detect(name, 'SREBF')) |>
  ggplot(aes(z_score, name)) +
  geom_col() +
  labs(title = '293 TG')

tg_293_xcore_chipatlas |>
  filter(str_detect(name, 'ATF4')) |>
  ggplot(aes(z_score, name)) +
  geom_col() +
  labs(title = '293 TG')

tg_hesc_xcore |>
  filter(str_detect(name, 'ATF')) |>
  ggplot(aes(z_score, name)) +
  geom_col()

# TRRUST -------
trrust <- read_tsv('00_util_scripts/ref/trrust_rawdata.human.tsv',
                   col_names = c('TF', 'target', 'type', 'ref'))

upr_target <- read_tsv('mission/FPP/thapsigargin_upr/upr_3tf_targets.tsv')

atf4_acti <- trrust |>
  filter(TF == 'ATF4', type == 'Activation') |>
  distinct(target) |>
  pull(target)

atf4_agi <- c('ATG7','HMOX1','SLC7A11','PSAT1','SESN2')

atf6_trrust <- trrust |>
  filter(TF == 'ATF6', type == 'Activation') |>
  pull(target) |>
  c(upr_target$ATF6f) |>
  unique()

xbp1s_trrust <- trrust |>
  filter(TF == 'XBP1', type == 'Activation') |>
  pull(target) |>
  c(upr_target$XBP1s) |>
  unique()

kera_infct_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera_v3_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3h_SA-PBS_deg.csv')

kera_v3l_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3l_SA-PBS_deg.csv')

hesc_tg_deg <-
  read_csv('mission/FPP/thapsigargin_upr/hesc_tg_vs_dmso_deg.csv') |>
  mutate(gene = Gene.Name, avg_log2FC = logFC, p_val_adj = FDR, .keep = 'none')

tg_293_deg <-
  read_csv('mission/FPP/thapsigargin_upr/293hek_tg_vs_dmso_deg.csv') |>
  mutate(gene = Name, avg_log2FC = logFC, p_val_adj = FDR, .keep = 'none')

uvb_skin_deg <- read_csv('mission/FPP/uvb/human.skin.bulk.6h.up.deg.csv') |>
  mutate(gene = Symbol, avg_log2FC = logFC, p_val_adj = FDR, .keep = 'none')

upr_sum <- list('SA KC' = kera_infct_deg, 'SA V3h KC' = kera_v3_deg,
     'SA V3l KC' = kera_v3l_deg, 'UVB skin' = uvb_skin_deg,
     'TG 293HEK' = tg_293_deg, 'TG hESC' = hesc_tg_deg) |>
  list_rbind(names_to = 'Condition') |>
  mutate(gene = str_to_upper(gene))

upr_sum |>
  write_csv('mission/FPP/thapsigargin_upr/upr6condition_tpm.csv')

upr_sum <-
  read_csv('mission/FPP/thapsigargin_upr/upr6condition_tpm.csv')

upr_atf4_sum <- upr_sum |>
  filter(gene %in% c(atf4_acti,atf4_agi))

upr_atf4_sum |>
  filter(!str_detect(Condition, 'V3')) |>
  write_source_csv('ATF4_UPR_condition_logfc')

## annotate sub-pathway ------------
atf4_go_anno <- upr_atf4_sum$gene |>
  unique() |>
  groupGO('org.Hs.eg.db', keyType = 'SYMBOL', ont = 'MF', readable = T, level = 4)

atf4_go_anno@result |>
  as_tibble() |>
  filter(Count > 0, str_detect(Description, 'popto|utophag|oxy|mino')) 

atf4_subpaths <- atf4_go_anno@result |>
  as_tibble() |>
  filter(Count > 0, ID %in% c('GO:0097190','GO:0006914','GO:0006520')) |>
  separate_longer_delim(geneID, delim = '/') |>
  mutate(pathway = Description, gene = geneID, .keep = 'none')

atf4_manual_path <-
  read_tsv('mission/FPP/thapsigargin_upr/atf4_downstream_pathways.tsv')

atf4_manual_path <- atf4_manual_path |>
  pivot_longer(-1, values_drop_na = T, names_to = 'Pathway')

atf4_kc_order <- upr_atf4_sum |>
  filter(Condition == 'SA KC') |>
  mutate(order = fct_reorder(gene, avg_log2FC), gene, .keep = 'none') 

g1 <- upr_atf4_sum |>
  left_join(atf4_kc_order) |>
  filter(!is.na(order), str_detect(Condition, 'V3', negate = T),
         str_detect(gene, 'FGF|NDC', negate = T)) |>
  ggplot(aes(Condition, order, fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(upr_atf4_sum$avg_log2FC)) +
  theme_pubr(base_size = 6, base_family = 'ArialMT', x.text.angle = 90,
             legend = 'right') +
  labs(title = 'ATF4-activated genes in UPR inducing condition', y = 'Gene') +
  theme_jpub

g2 <- atf4_manual_path |>
  left_join(atf4_kc_order) |>
  filter(!is.na(order)) |>
  ggplot(aes(Pathway, order, fill = Pathway)) +
  geom_tile() +
  theme_pubr(base_size = 6, base_family = 'ArialMT', legend = 'right') +
  labs(y = '') +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
  theme_jpub

g1 + g2 + patchwork::plot_layout(guides = 'collect')

publish_pdf('atf4_downstream_2heatmap.pdf', width = 80, height = 70)

upr_atf6_sum <- upr_sum |>
  filter(gene %in% atf6_trrust)

upr_atf6_sum |>
  filter(Condition == 'SA V3h KC') |>
  mutate(order = fct_reorder(gene, avg_log2FC), gene, .keep = 'none') |>
  right_join(upr_atf6_sum) |>
  filter(!is.na(order)) |>
  ggplot(aes(Condition, order, fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(upr_atf6_sum$avg_log2FC)) +
  theme_pubr() +
  labs(title = 'ATF6-activated genes in UPR inducing condition', y = 'Gene')

upr_xbp1_sum <- upr_sum |>
  filter(gene %in% xbp1s_trrust)

upr_xbp1_sum |>
  filter(Condition == 'SA V3h KC') |>
  mutate(order = fct_reorder(gene, avg_log2FC), gene, .keep = 'none') |>
  right_join(upr_xbp1_sum) |>
  filter(!is.na(order)) |>
  ggplot(aes(Condition, order, fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdBu',
                       values = pretty_distiller(upr_xbp1_sum$avg_log2FC)) +
  theme_pubr() +
  labs(title = 'XBP1s-activated genes in UPR inducing condition', y = 'Gene')

## export tpm ---------
atf4_selected <-
c('SLC7A11','DDR2','SIGMAR1','PSAT1','ASNS','MAP1LC3B','ATG5','VEGFA','ATG7')

hek293 |>
  pivot_longer(-1) |>
  calc_tpm(sample = name, abundance = value) |>
  filter(Name %in% atf4_selected) |>
  mutate(sample = name, gene = Name, log_tpm = log1p(tpm), tpm, .keep = 'none') |>
  write_source_csv('TG_hek293_atf4_downstream_tpm')

hesc_tdb |>
  calc_tpm(sample = name, abundance = value) |>
  filter(Gene.Name %in% atf4_selected) |>
  mutate(sample = str_c(group, suffix), gene = Gene.Name, 
         log_tpm = log1p(tpm), tpm, .keep = 'none') |>
  write_source_csv('TG_hESC_atf4_downstream_tpm')

tdb.sk |>
  filter(Symbol %in% atf4_selected) |>
  mutate(sample = str_extract(title, '\\d+_hour.+'), gene = Symbol, 
         log_tpm = log1p(tpm), tpm, .keep = 'none') |>
  write_source_csv('UVB_skin_atf4_downstream_tpm')
